Gamma cameras known in the art of Nuclear Medicine (NM) imaging produce tomographic images that are indicative of physiological activity. Such cameras receive radiation that is emitted by radioisotope markers or tracers, which are introduced into the body of a subject and are taken up by an organ of the body in proportion to the physiological activity of interest. The radiation emitted is generally received by a scintillator/detector system, which produces electrical signals responsive to photons of the radiation incident thereon. The signals are processed and back-projected, using computerized tomography methods known in the art, to produce a three-dimensional image indicative of localized activity within the organ.
Positron emission tomography (PET) is a system of tomographic nuclear imaging which is well known. In general, this system is based on utilization of radio-isotopes which, during a decay, emit two photons in directly opposite directions. A ring, or rings of discrete detectors which surround a subject into whose body such isotopes have been introduced, detects the occurrence of such a decay by detecting two coincident gamma rays impinging at two detectors, where the events have an energy associated with the decay.
Based on this coincidence detection, the position of the decay (i.e., the presence of the radio-isotope which decayed and caused the coincident detection) is known to be along the line joining the two positions at which the coincident impingements were detected.
The exact origin of the event is not known and the calculation of an image of the distribution of the radio-isotope is, in the prior art, based on a probabilistic smearing of the event into sinograms associated with the line on which the event is known to have occurred. This smeared probability is forward projected to form (together with other detected events) tomographic views at each of a plurality of slice positions, generally coincident with the rings of detectors. These views are used to generate tomographic images of slices corresponding generally to the positions of the rings.
One of the major problems with this reconstruction system is the “blurring” of events in the direction perpendicular to the plane of the rings. This blurring, if not reduced, results in an image which is unsuitable, for diagnostic purposes.
One method of reducing the effects of axial blurring is just to reduce the “acceptance angle” for events. If events cause coincidence detection in widely spaced rings, they naturally cause greater axial blurring. By reducing the acceptance angle, i.e., the angle of the line connecting the detection points with a ring, amount of axial blurring may be reduced, at the price of rejecting valid detected events.
A second method, described in U.S. Pat. No. 5,331,553, to Muehllehner et al., axially rebins the events based on a deblurring function. This rebinding reduces the axial blurring, however, it increases the image noise. Substantial blurring remains and artifacts are generally generated. The rebinning may be performed on the views or may be performed on the three dimensional tomographic image.
Another method which has been mooted is to use an expectation maximization (EM) algorithm 3-D reconstruction. Such an algorithm is described, for example, in “Maximum Likelihood Reconstruction for Emission Tomography” IEEE Trans Med Imag MI-1; pages 113–122 (1982). This algorithm performs iterative expectation maximization operations on an equation relating the activity in the voxels with the projection data. During each iteration, all of the projection data is taken into account. However, this method, while theoretically possible, requires very large amounts of computation in order to reach satisfactory results, so that it has not been implemented commercially.
In order to more clearly understand the operation of the invention, it is useful to review the prior art EM methods.
In the EM method, a body having a variable radiation emission density is considered to be contained in a discrete cube or other discrete region and the emission density of each voxel, v=(x,y,z) in the region is defined as λ(v). The radiation emitted by the body is detected by detectors surrounding the body. If two photons are simultaneously detected by the two detectors (which indicates that the event may be caused by a positron interaction) and a Line Of Flight (LOF) between the two detection coordinates intersects the region we consider this as a coincidence acquisition (more simply as an “event”). Coincidence acquisition reconstruction algorithms try to determine the unknown emission density distribution in the region given a list of the LOF of detected events.
The classic PET scanner is built of rows of detector rings whereby each detector is a discrete unit. All the events detected simultaneously by the same two detectors di and dj are collected in a single bin bij. Thus each bin b defines a single LOF and every event detected in the bin is assumed to have originated along this LOF.
Let v=1, . . . , V represent voxels of the field of view and let independent variables x(v) with unknown mean values of emission density λ(v) represent the number of unobserved emissions in each of the V voxels. Suppose further that if an emission occurs in voxel v it has a probability of p(v,b) of being detected in bin b, then p(v,b) defines a transition matrix (likelihood matrix) which is known. Based on the number of y=y(b) events detected in each bin it is desired to estimate the total number of the unknown distribution of events λ=λ(v), v=1, . . . , V. For each λ, the observed data has the conditional probability or likelihood:
                    P        (                              y            ⁢                                        λ              )                                =                                    ∏                              1                ≤                b                ≤                B                                                                                  ⁢                                                  ⁢                                          ⅇ                                                      -                    μ                                    ⁢                                                                          ⁢                  b                                            ⁢                                                                    μ                    ⁡                                          (                      b                      )                                                                            y                    ⁡                                          (                      b                      )                                                                                                            y                    ⁡                                          (                      b                      )                                                        !                                                                                        (        1        )            where μ(b) are the mean values of the (observed) Poisson variables y(b), that is:
                              μ          ⁡                      (            b            )                          =                              ∑                          1              ≤              v              ≤              V                                                                      ⁢                                          ⁢                                    λ              ⁡                              (                v                )                                      ⁢                                          p                ⁡                                  (                                      v                    ,                    b                                    )                                            .                                                          (        2        )            
The maximum likelihood estimate of λ is:{overscore (λ)}=arg maxλP(y/λ)  (3)
The optimality condition of the above equation, based on the log-likelihood, is:
                                          λ            ⁡                          (              v              )                                =                                    1                              P                ⁡                                  (                  v                  )                                                      ⁢                                          ∑                                  b                  =                  1                                B                            ⁢                                                          ⁢                                                                    y                    ⁡                                          (                      b                      )                                                        ⁢                                      λ                    ⁡                                          (                      v                      )                                                        ⁢                                      p                    ⁡                                          (                                              v                        ,                        b                                            )                                                                                                            ∑                                                                  v                        ′                                            =                      1                                        V                                    ⁢                                                                          ⁢                                                            λ                      ⁡                                              (                                                  v                          ′                                                )                                                              ⁢                                          p                      ⁡                                              (                                                                              v                            ′                                                    ,                          b                                                )                                                                                                                                ,                  v          =          1                ,        …        ⁢                                  ,        V                            (        4        )            where P(v) is the probability to detect an event emitted from voxel v:
                              P          ⁡                      (            v            )                          =                              ∑                          b              =              1                        B                    ⁢                                    p              ⁡                              (                                  v                  ,                  b                                )                                      .                                              (        5        )            
The EM algorithm can be considered as a fixed point iterative algorithm based on (4):
                                                        λ                              n                ⁢                                                                  ⁢                e                ⁢                                                                  ⁢                w                                      ⁡                          (              v              )                                =                                    1                              P                ⁡                                  (                  v                  )                                                      ⁢                                          ∑                                  b                  =                  1                                B                            ⁢                                                                    y                    ⁡                                          (                      b                      )                                                        ⁢                                                            λ                                              o                        ⁢                                                                                                  ⁢                        l                        ⁢                                                                                                  ⁢                        d                                                              ⁡                                          (                      v                      )                                                        ⁢                                      p                    ⁡                                          (                                              v                        ,                        b                                            )                                                                                                            ∑                                                                  v                        ′                                            =                      1                                        V                                    ⁢                                                            λ                                              o                        ⁢                                                                                                  ⁢                        l                        ⁢                                                                                                  ⁢                        d                                                              ⁢                                                                                  ⁢                                          (                                              v                        ′                                            )                                        ⁢                                          p                      ⁡                                              (                                                                              v                            ′                                                    ,                          b                                                )                                                                                                                                ,                  v          =          1                ,        …        ⁢                                  ,        V                            (        6        )            
In 1994, a paper entitled “Accelerated Image Reconstruction Using Ordered Subsets of Projection Data” IEEE Trans Med Imag vol. 13, no. 4, (1994) pp. 601–609 reported the use of an EM algorithm using ordered sub-sets which it called OSEM. The OSEM algorithm also performs an EM algorithm on the equation relating the voxel activity with the projection data. However, in this method sub-sets of the projection data, rather than the full set of projection data, are taken into account in each iteration with a different sub-set being taken into account for each iteration. If the projection data is divided into N sub-sets (which together form the complete projection data set) and each iteration is performed using only one of the sub-sets, then, if N iterations of this type are performed, the overall results will be comparable to those achieved when N iterations are performed taking into account all of the projection data.
This surprising result allows the practical use of EM reconstruction algorithms (in the OSEM form) in PET systems using rings of detectors. A description of the OSEM method as in the above referenced paper by Hudson, et al., follows:
Let yi be the number of photon emissions recorded in the ith projection bin and let Yθ be the set of parallel projections {yi, yi+1, . . . } that view the object at angle θ, orthogonal to the tomographic axis. The projection data are grouped into n subsets S1, S2, . . . ,Sn. If there are a total of P projection bins, the elements of each subset, in no particular order, are:S1={y1,y2, . . . yP/n}, S2={y(P/n)+1,y(P/n)+2, . . . y2P/n}, . . . ,Sn={y(n−1)(P/n)+1,y(n−1)(P/n)+2, . . . yp}.  (7)
The subsets normally consist of projection views separated by some fixed angle about the object. For example, each subset might consist of two sets of parallel projection views spaced 90° apart, e.g., S1={(Y0, Yπ/2}, S2={(Yπ/4,Y3π/4}, and so on.
In the OSEM algorithm, the log-likelihood objective function for each of the subsets is increased with each iteration, using the results of the previous iteration as the starting point. Therefore, the EM iterations become:
                                                        λ                              k                +                1                                      ⁡                          (              v              )                                =                                    1                              P                ⁡                                  (                  v                  )                                                      ⁢                                          ∑                                  b                  ⁢                                                                          ⁢                  ε                  ⁢                                                                          ⁢                                      S                    k                                                                                                                ⁢                                                                    y                    ⁡                                          (                      b                      )                                                        ⁢                                                            λ                      k                                        ⁡                                          (                      v                      )                                                        ⁢                                      p                    ⁡                                          (                                              v                        ,                        b                                            )                                                                                                            ∑                                                                  v                        ′                                            =                      1                                        V                                    ⁢                                                            λ                      k                                        ⁢                                                                                  ⁢                                          (                                              v                        ′                                            )                                        ⁢                                          p                      ⁡                                              (                                                                              v                            ′                                                    ,                          b                                                )                                                                                                                                ,                  v          =          1                ,        …        ⁢                                  ,                  V          .                                    (        8        )            where λk is the estimated number of emissions from v after the introducing the kth subset of projections.
The EM procedure is repeated until all n subsets have been exhausted. Such a cycle is considered a single iteration of OSEM. The cycle can be repeated iteratively until a satisfactory reconstruction is obtained.
However, ring type PET systems are not of general utility. Their use of separate detectors for each detection pixel results in a system which can have high sensitivity, but not high resolution. For this reason and, to a lesser extent, because of its geometry, such a system is useful only for PET and cannot be used for other NM applications such as for acquiring planar images and for SPECT.
Neither the EM nor OSEM algorithms are easily applied to systems of planar detectors. Such application would be very desirable since this would allow for the use of standard rotating dual head gamma cameras for PET as well as for SPECT. However, this is not practical. If an attempt is made to use a pair of opposed planar detectors (of the type normally used for planar or SPECT imaging), rotating about the subject, to acquire data for constructing PET images using EM (or OSEM) techniques, the amount of data which is acquired is reduced, as compared to the ring type PET system, making the planar system impractical. For example, if two detectors having dimensions of 540×400 mm and a rotation radius of 350 mm are used, with a bin resolution of 2.5×2.5 mm, there are about 3×108 bins. If the system is capable of an acquisition rate (for coincident photon events) of 1000 events per second, in a typical study of 30 minutes approximately 1.8 million events can be collected. The number of bins is more than two orders of magnitude larger than the number of events, which means that most bins will remain empty and almost no bins will have multiple events. Thus, for such a device, the normal binning procedure according to a fixed set of a very large number of bins is impractical and results in excessive calculation. Also, with such sparse data the results can be expected to be noisy.
It is understood that for the classical ring type PET imager, the number of bins is much smaller and the absorption of the photons is more efficient (since the crystal is thicker) such that, OSEM reconstruction becomes practical. The prior art does not teach any way to combine the advantages of the OSEM system with the higher intrinsic resolution of planar detectors.